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A three-dimensional Eulerian analysis has been developed for modeling droplet 
impingement on lifting bodes. The Eulerian model solves the conservation equations of mass 
and momentum to obtain the droplet flow field properties on the same mesh used in CFD 
simulations. For complex configurations such as a full rotorcraft, the Eulerian approach is 
more efficient because the Lagrangian approach would require a significant amount of 
seeding for accurate estimates of collection efficiency. Simulations are done for various 
benchmark cases such as NACA0012 airfoil, MS317 airfoil and oscillating SC2110 airfoil to 
illustrate its use. The present results are compared with results from the Lagrangian 
approach used in an industry standard analysis called LEWICE. 


I. Introduction 

I n-flight icing of a helicopter is an interdisciplinary problem and an active area of concern. Icing affects the 
availability, safety and survivability of the rotorcraft due to degradation in performance caused by loss of lift, drag 
increase, and decrease in stall angle. A first step in modeling ice accretion is to estimate the rate at which water 
droplets collect on the solid surface. Advances in understanding and estimating water droplet impingement on 
aerodynamic bodies have been made through experimental studies and through the development of analytical and 
empirical tools over the past few decades. 1-3 

There have been two primary approaches for the prediction of surface droplet impingement distributions - 
Lagrangian and Eulerian methods. Da Silveira et al 4 have conducted an evaluation of these methods and found both 
methods to be equally effective. LEWICE or LEWICE3D are representative examples of industry- standard icing 
programs that use a Lagrangian approach to compute droplet trajectories through the air, and have been shown to be 
highly effective. 5,6 In Lagrangian approaches, computational cost is reduced by performing the simulation of ice 
accretion only at a few selected strips in the configuration, as opposed to the full 3D simulation where collection 
efficiency is computed over the entire surface. 

In the Eulerian approach, (e.g. FENSAP-ICE. 7,8 ) the conservation of mass and momentum of the droplets are 
computed simultaneously with the flow field solution, by solving two additional governing equations for the volume 
fraction of water and the particle velocities. These equations are solved on the same CFD mesh. The mean flow may 
be unsteady, and the solid surfaces may be in relative motion. Most Lagrangian approaches, on the other hand, 
assume or require the flow field to be steady. For this reason, an Eulerian method is more attractive for modeling 
rotorcraft icing phenomena. 

The present researchers, in collaboration with NASA and industries, have been working on the development of 
methods for modeling rotorcraft icing phenomena. In the past, this effort was based on Lagrangian approaches. 9 ' 14 
The primary objective of the present study is to couple an existing CFD analysis with an Eulerian droplet solver and 
to verify the accuracy and the computational efficiency of this method for modeling droplet impingement. In order 
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to validate the present approach using Eulerian method, simulations are performed for various benchmark cases. The 
present computational results are compared with results from LEWICE simulations and test data. 

II. Icing model formulation 

Figure 1 shows the basic elements of the ice accretion simulation model. The process starts with grid generation 
and CFD flow analysis for a clean baseline configuration. The droplet solver reads the flow field data and computes 
the local collection efficiency (p) on the surface. This information is fed into LEWICE which subsequently 
calculates the resulting ice shape that evolves over a period of time. The grid generator is next invoked to generate a 
new volume grid around the iced configuration, for use in the CFD solver for an updated flow filed. These modules 
are coupled to each other using a PYTHON script, and exchange the required data using industry -standard flow filed 
and grid format (PLOT3D). In the present study the focus is mainly on the development and validation of the 
Eulerian droplet solver. 


Clean Geometry 


CFD Code 
(Flow solution) 


(Build grid for iced rotor) 


Droplet Solver 
(Collection efficiency) 


LEWICE 
[ Ice Shape) 


Figure 1. Overview of the Ice Accretion Analysis 


A. CFD solver 

A 3-D Reynolds -Averaged compressible Navier-Stokes solver called GENCAS (Generic Numerical 
Compressible Airflow Solver), developed by Min 15, 16 is used in this study to model the flow field. This solver may 
be used to solve the RANS equations on 2D or 3D structured multi -block grids. Roe’s FDS and AUSMPW+ upwind 
schemes are available for inviscid flux computation. First or second order implicit LUSGS with Newton sub- 
iteration, or 2 nd / 4 th order explicit Runge-Kutta schemes are available for time marching. If a higher order spatial 
accuracy is desired, 3rd order MUSCL, 5 th order WENO or a 7 th order WENO cell interface reconstruction methods 
may be selected. Available turbulence models include one equation Spalart-Allmaras (SA) and SA-DES models, and 
two equation Wilcox’s k-cd, standard k-s, Menter’s k-cd/k-s BSL, Menter’s k-co SST, KES, and HRKES models. For 
a detailed description of the numerical formulation of GENCAS, the reader is referred to the papers written by Min 
etal.. 15 ’ 16 

B. Droplet solver 

In order to compute the droplet flowfield properties at the same nodes of the discrete domain where the flow 
variables of air are known, an Eulerian approach is used in the present study. In this method, the average water 
droplet properties within a control volume are solved instead of tracking individual particles. This physical approach 
has several advantages over the Lagrangian approach. These include improved quality of the solution, the ability to 
model unsteady flows over bodies in relative motion, and the automated treatment of shadow zones (no 
impingement) for probes or detector placing. 17 The interaction between the air particles and the droplets occurs 
through a drag force exerted by the mean flow on the particles. The presence of the droplet flow field is not felt by 
the mean flowfield solver, and the droplets are treated as a passive scalar field. When the air flow is steady, the CFD 
analysis may be computed a priori and used in the droplet solver. 

In the derivation of governing equations for air-droplet flows, the following assumptions are made 18 : 

2 

American Institute of Aeronautics and Astronautics 






■ The droplets have a spherical shape and do not undergo any deformation or breaking. 

■ There is no collision, or coalescence between droplets. 

■ There is no exchange of heat and mass between the droplets and the surrounding air. 

■ The effect of mean flow mixing effects on the droplet is neglected. 

■ Drag, gravity and buoyancy due to density differences are the only forces acting on the droplets. 

The first two assumptions are based on the fact that the size of icing droplets is 1-100 pm range and droplet flow 
is considered dilute with a volume fraction around 10‘ 6 . Although the gravity and buoyancy forces are three orders 
lower in magnitude than drag force in typical flight conditions, these forces are kept in the model because their 
effects could be significant in the simulation of de-icing fluid contamination by rain and snow during ground 
operation. 

The governing equations for the conservation of mass and momentum of the droplets are written as follows: 

^ + V-(« M ,.) = 0 (1) 

at 
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Here, u a is the non-dimensionalized velocity of air; p ,the density of water; p a ,the density of air; g h gravity 


vector; F r = U^I -jLg is the Froude number; ,the speed of air at freestream; L, the characteristic length 


(typically the airfoil chord length); K = pd 2 U ^ H%Lp , an inertia parameter; p ,the dynamic viscosity of air. 

The first term on the right-hand- side of the momentum equation accounts for the drag acting on the droplet or 
particle based on low-Reynolds number, or Stokes flow, behavior for spheres. 19 The droplets Reynolds number (Re d ) 
is defined based on the slip velocity between the air and droplet and the droplet diameter. The drag coefficient is 

C D = -^(l + 0.15Re/ 687 ) Re d <1000 

Re/ ’ (3) 
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Equations (1) and (2) are recast in finite volume form using divergence theorem. A first order upwind scheme is 
employed for computing the mass and momentum flux at the faces of the control volume. An implicit time marching 
algorithm is employed. Mean flow quantities are lagged by one time step compared to particle velocity and volume 
fraction. The resulting system of banded linear equations is solved using an approximate factorization scheme. 

The freestream values of droplet velocity and volume fraction are imposed as boundary conditions at the far 
field. Prescribing the correct boundary conditions for the droplets at the wall is not straight-forward. The droplet 
velocity cannot be simply set to zero on the walls. A switching boundary condition 20 is applied. Volume fraction and 
velocity of droplets are extrapolated from the computed flux entering the control volumes adjacent to the solid. A 
lower bound of volume fraction and zero velocity are imposed on flux exiting the flowfield, and collecting on the 
walls. 


a w =a w -l andui w = u i, w -i Incomingdropletfluxes 
tr w = f a and^ w =0 Dropletfluxexitingthewall 


A common way of comparing droplet impingement rate at various flight conditions is through the collection 
efficiency (p). This quantity characterizes the configuration's ability to capture incoming water and is defined as the 
local mass flux of water onto the airfoil surface normalized by the freestream liquid water content and the freestream 
velocity. 


. = _Wi_A 

(LWCpn |A| 


(5) 


where, A { is the local area normal; LWC , Liquid water content (Fig. 2). 

Computational prediction for large droplet case and found to show a considerably higher collection efficiency 
distribution and the peak value is greater than the measurement. 21 A plausible reason for this over prediction is 
droplet splashing and breakup. 21 In present study, the effect of droplet splashing is investigated by using a model 
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proposed in Ref. 22. Splashing causes a reduction in collection efficiency. The mass faction of water lost due to 
bouncing (N b ) is first computed. 
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Here, c> is surface tension between air and water. Finally the collection efficiency is computed as 

p'=p{\-N h ) 
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Figure 2. Definition of Collection Efficiency 

III. Results 

In this section, a number of calculations are presented to demonstrate the capabilities of the present Eulerian 
approach. Comparisons with industry- standard Lagrangian approaches found in LEWICE are also shown. 

A. NACA0012 airfoil 

Collection efficiency estimates have been done for NACA0012 airfoil, at three different angles of attack. The 
simulations are performed at a 0.31 free-stream Mach number with a constant droplet diameter of 20pm and an 
airfoil chord of 0.5334 m. The mean flow field is obtained from GENCAS. 15-16 In the CFD simulation, Roe scheme 
with a 3rd order MUSCL reconstruction is used for flux calculations. A first order implicit LU-SGS scheme is used 
for marching in time. Spalart-Allmaras (SA) is used as the turbulence model. 

Figure 3 shows a comparison of the collection efficiency from the present Eulerian simulations with the 
LEWICE Lagrangian results. 20 In general, the present results are in good agreement with LEWICE, providing 
confidence in the present method. It is found that the deviation between the two approaches grows with increased 
angles of attack. Similar discrepancies have been reported by Kinzel et al. 20 and Beaugendre et al. 23 in their 
comparisons between FENCAP-ICE and LEWICE. For the a =4 deg. the results from LEWICE are obtained at 
corrected angle of attack (3.5 deg.). 


Table 1. Test Conditions for MS317 Airfoil 


Parameter 

Value 

Chord (m) 

0.914 

Uoo (m/sec) 

78.66 

Re (Million) 

4.83 

AOA (Degree) 

0/8 

MVD (pm) 

11.5/21.0/92.0 
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B. MS317 airfoil 

Next, the collection efficiency simulations are reported for the MS317 airfoil. This configuration was chosen 
because of the availability of collection efficiency and pressure distributions data at various mean flow conditions, 
collected over 1997 and 1999. 24 GENCAS is used to obtain flowfield data. In the CFD simulation, Roe scheme with 
a 3rd order MUSCL reconstruction is used for flux calculations. 1 st order implicit LUSGS scheme is used for 
marching in time. Spalart-Allmaras (SA) is used as a turbulence model. The predicted pressure distributions are 
compared with experimental data in Fig. 4. Predicted pressure distributions at the bottom surface are in good 
agreement with experiment. Some differences between the computed and measure pressure distributions are 
observed near the trailing edge, but this is expected to play on a minor role in the collection efficiency near the 
leading edge. 

The effect of median volumetric diameter (MVD) on collection efficiency is investigated. The icing test 
conditions are given in table 1. The effect of first cell distance normalized by chord length was also investigated 
because the droplet solver updated the values at boundary by using the values of first inner cell. It is found that the 
collection efficiency is relatively insensitive to the normal height of the first row of cells over the wall. It is expected 
that the deviation in the flowfield between present simulation and the test data would only have a negligible effect 
on the collection efficiency distribution around the leading edge. In the experiment, collection efficiency was 
measured for 0 and 8 degree of angle of attack and MVDs of 1 1.5, 21, and 92pm. 

Figure 5, 6 present the comparison of local impingement efficiency distributions between present prediction and 
measurement according to different value of MVD at 0° and 8°. The x-axis (surface distance) is normalized by 
airfoil chord length. The positive values correspond to the lower surface of the airfoil. The peak value of collection 
efficiency increases with MVD size. For an angle of attack of 0°, the laser system shows higher impingement 
efficiency values near the region of maximum impingement efficiency. In Ref. 24, the reason for this discrepancy is 
explained. It was attributed to a small level of dye penetration into the blotter. In the present simulation, the 
impingement limits are under-predicted except for the 92 pm case for which the predicted collection efficiency is 
considerably higher and the peak value is greater than the measurement. A similar over -prediction is seen in the 
results from FEWICE in Ref. 24. Possible reasons for these large differences between simulation and experiment 
was investigated in Ref. 24. One of the cited reasons was the errors associated with measuring MVD for the 92-94 
pm cases. Another plausible reason is droplet splashing and breakup. 

Additional studies were performed for this test condition and it was found that droplet splashing and breakup 
occurs near the airfoil leading edge region. For the high angle of attack case, the location of peak value of collection 
efficiency was shifted downstream on the lower surface of the airfoil. Simulation results are shifted to the left with 
respect to the experimental data, if the angle of attack is not corrected for wall effects. 

The effect of first cell distance on collection efficiency is investigated in Fig 5 -(b). Marginal difference in 
collection efficiency is observed. The effect of droplet splashing is investigated in Figures 5 -(c) and 6-(c). An 
improvement in the prediction is seen when the collection efficiency is modified to account for splashing. 





(a) a =0 deg. (b) a =4 deg. (c) a =8 deg. 

Figure 3. Comparison of Predicted versus LEWICE Collection Efficiencies for a NACA0012 Airfoil (Wall 

Corrections have not been used in the simulations) 
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Figure 4. Comparison of Pressure Distribution for MS317 Airfoil 




Surface Distance from Highlight Surface Distance from Highlight 

(a) MVD=11.5 (b) MVD=21 



Surface Distance from Highlight 


(c) MVD=92 


Figure 5. Comparison of Collection Efficiency for MS317 airfoil at Zero Angle of Attack 
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(a) MVD=11.5 (b) MVD=21 



(c) MVD=92 


Figure 6. Comparison of Collection Efficiency for MS317 airfoil at 8 Degrees Angle of Attack 
C. Oscillating SC2110 airfoil 

Collection efficiency calculations have been performed for an oscillating SC2110 airfoil and comparisons with 
LEWICE have been made. The airfoil has a chord length of 0.381m, and operates at a freestream Mach number of 
0.4208. Unsteady flowfield data for each angle of attack were obtained using a version of OVERFLOW. 

The simulations employ a nominal MVD size of 22 pm. The collection efficiency is computed for -1, -0.75, 0.15, 
5, 8.53 and 11 degrees of angle of attack. Comparisons of collection efficiency between the present simulation and 
LEWICE for oscillating SC2110 airfoil are presented in Fig. 7 at several angles of attack. The present Eulerian 
approach shows a spatial distribution of collection efficiency similar to LEWICE. The peak values from the two 
approaches are in reasonable agreement. It is found that the present Eulerian simulation shows a wider surface 
region with significant collection of water droplets compared to the Lagrangian simulation. 
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Collection Efficiency (p) Collection Efficiency (p) Collection Efficiency (p) 



Normalized Arc Length (s/c) 

(a) a =-l deg. 





Normalized Arc Length (s/c) 


Normalized Arc Length (s/c) 


(c) a =0.15 deg. 


(d) a — 5 deg. 




Normalized Arc Length (s/c) Normalized Arc Length (s/c) 

(e) a =8.53 deg. (f) a = 11 deg. 

Figure 7. Comparison of Collection Efficiency for an Oscillating SC2110 Airfoil 
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Table 2. Test Conditions for NACA64A008 Swept Tail 


Parameter 

Value 

Uoo (m/sec) 

78.66 

Re (Million) 

5.03 

AOA (Degree) 

0/6 

MVD (pm) 

11.5/21.0 


D. NACA64A008 Swept Tail 

In an effort to assess the suitability of the present approach for 3-D configurations, collection efficiency 
simulations have been done for a swept tail made of NACA64A008 sections. This configuration was chosen because 
of the availability of collection and pressure distributions data at various mean flow conditions, collected over 1997 
and 1999. 24 GENCAS is used to obtain flowfield data. The predicted pressure distributions are compared with 
experimental data in Fig. 7 and are in good agreement with experiment. 

The icing test conditions are given in Table 2. Figure 9 and 10 present the comparison of local impingement 
efficiency distributions between present prediction and measurement according to different value of MVD at 0° and 
6°. The x-axis (surface distance) is normalized by airfoil chord length. The positive values correspond to the lower 
surface of the tail section. The peak value of collection efficiency is found to increase with MVD size. For an angle 
of attack of 0°, the peak values of collection efficiency are under-predicted. For the high angle of attack case, the 
location of peak value of collection efficiency was shifted downstream on the lower surface of the airfoil. 
Simulation results are shifted to the left with respect to the experimental data. 

IV. Conclusions 

A 3D Eulerian based stand-alone solver has been developed and tested for various benchmark cases. The present 
Eulerian based solver has been shown to successfully predict collection efficiencies on two-dimensional and quasi 
three dimensional airfoils. The present approach is also in reasonable agreement to a well -validated Fagrangian code 
(FEWICE). 
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a= 0 deg, M=0.23, Re = 5.03 million 




(a) a =0 deg. 

Figure 8. Comparison of Pressure Distribution for NACA64A008 Swept Tail Section 
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(a) MVD=11.5 (b) MVD=21 

Figure 9. Comparison of Collection efficiency for NACA64A008 Swept Tail Section at Zero Degrees 

Angle of Attack 



Figure 10. Comparison of Collection Efficiency for NACA64A008 Swept Tail at 6 Degrees Angle of Attack 

(MVD=21) 
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